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Abstract 

We conduct an investigation into the dispersive post-shock oscillations in the entropic lattice- 
Boltzmann method (ELBM). To this end we use a root finding algorithm to implement the ELBM which 
displays fast cubic convergence and guaranties the proper sign of dissipation. The resulting simulation 
on the one- dimensional shock tube shows no benefit in terms of regularization from using the ELBM over 
the standard LBGK method. We also conduct an experiment investigating of the LBGK method using 
median filtering at a single point per time step. Here we observe that significant regularization can be 
achieved. 

Keywords: Fluid dynamics, lattice Boltzmann, entropy balance, dispersive oscillations, numerical test, 
shock tube 

AMS subject classifications. 65N12, 76M28, 74Q10, 74J40 

1 Introduction 

The Lattice Boltzmann methods in their original form (see [5J I12| ) do not guarantee the proper entropy 
production and may violate the Second Law. The proper entropy balance remains up to now a challenging 
problem for many lattice Boltzmann models [14] . 

The Entropic lattice Boltzmann method (ELBM) was invented first in 1998 as a tool for construction of 
single relaxation time lattice Boltzman models which respect the TJ-theorem [9]. For this purpose, instead 
of the mirror image with local equilibrium as reflection center, the entropic involution was proposed, which 
preserves the entropy value. Later, we call it the Karlin-Succi involution [7]. In 2000, it was reported 
that exact implementation of the Karlin-Succi involution (which keeps the entropy balance) significantly 
regularizes the post-shock dispersive oscillations [T|. This regularization seems very surprising, because the 
entropic lattice BGK (ELBGK) model gives a second-order approximation to the Navicr-Stokes equation 
(different proofs of that degree of approximation were given in |12j and [1] ) , and due to the Godunov theorem 
[5] linear second-order finite difference methods have to be non monotonic. 

Moreover, Lax |10j and Levermore with Liu [llj . demonstrated that these dispersive oscillations are 
unavoidable in classical numerical methods. Schemes with precise control of entropy production, studied 
by Tadmor with Zhong [13], also demonstrated post-shock oscillations. Of course, there remains some gap 
between methods with proven existence of dispersive oscillations, and ELBM. However, recently, the existence 
of oscillations in the vicinity of the shock, at small values of viscosity for ELBM, was reported for Burgers' 
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equation [3]. In a recent paper [5] post shock oscillations of ELBGK were reported too, and no difference 
was found between ELBGK and LBGK in that regard. 

Nevertheless, absence of dispersive oscillations for ELBGK was reported many times since 2000. In this 
paper we answer the question: does the precise control of entropy production by ELBGK smooth the post- 
shock oscillation? The answer is negative. The exact implementation of the entropic involution does not 
smooth the dispersive oscillation (similarly, the exact control of entropy production does not smooth the 
post shock oscillation in finite difference methods [13]). Hence, the smoothing effect is caused by numerical 
imprecision in calculations of entropic involution, i.e. in solution of the following transcendental equation 
with respect to a (a^O): 

S(f + a(f*-f)) = S(f), (1.1) 

where S is entropy / is a current distribution, and /* is the corresponding equilibrium. 

In the first part of this paper we discuss a different numerical implementation of the ELBGK and conduct 
an investigation into exactly what stabilization properties it exhibits. 

The choice of the method for solution of (|1.1|) should be very precise, and in Section @] we describe a 
cubically converging root finding algorithm. It is not sufficient to have high precision when we have average 
deviation of the current distribution / from the associated equilibrium /*. For example, for solutions with 
shocks, it is usual for the distribution of this deviation to far from being exponential [5], and there appear 
points with deviation of several orders higher than the average. Moreover, it is sufficient to smooth a solution 
at one point only. We demonstrate this in the second part of the paper. We select the lattice site with most 
noncquilibrium / and regularize the field of noncquilibrium entropy at this point with 3-point median filter 
[5]. As a result, the dispersive oscillations drastically decrease. 



2 Lattice Boltzmann methods 

The Lattice Boltzmann method arises as a discretization of Boltzmann's kinetic transport equation 

^+v-V/ = Q c (/). (2.2) 

The population function / describes the distribution of single particles in the system and the collision integral 
Q c their interaction. Altogether (|2.2| describes the behaviour of the system at the microscopic level. By 
selecting a finite number of velocities v,;, (i = 1, . . . ,n) we create discrete approximation of the kinetic equation 
in velocity space. An appropriate choice of the velocities and time step discretizes space. For a time step 
of St= 1 the lattice can be created by unsealed space shifts of the velocities, and we get the fully discrete 
lattice Botzmann gas: 

fi(x + Vi,t+l) = fi(x,t) + Qi (2.3) 

where the proper transition from continuous collision integral Q c (f) to its fully discrete form {Qi} is assumed. 
The simplest and the most common choice for the discrete collision integral Qi is the Bhatnagar-Gross-Krook 
operator with over-relaxation 

Qi = a0(f*-fi). (2.4) 

For the standard LBGK method a = 2 and /3s[0,l] (usually /?£ [1/2,1]) is the over-relaxation coefficient 
used to control viscosity. For (5=1/2 the collision operator returns the local equilibrium f* and (3=1 (the 
mirror reflection) returns the collision for a liquid at the zero viscosity limit. For a viscosity v the parameter 
(3 is chosen by (3 = 8t/(2v + 8t). It should be noted that a collision integral such as (|2.4[) demands prior 
knowledge of a local equilibrium state for the given lattice. 
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A variation on the LBGK is the ELBGK pQ. In this case a is varied to ensure a constant entropy 
condition according to the discrete if-theorem. In general the entropy function is based upon the lattice and 
cannot always be found explicitly. However in the case of the simple one dimensional lattice with velocities 
v=(— c,0,c) and corresponding populations f=(/_,/o,/+) an explicit Boltzmann style entropy function is 
known [5]: 



The trivial root a = returns the entropy value of the original populations. ELBGK then finds the non- 
trivial a such that (|2.6[) holds. This version of BGK collision one calls entropic BGK (or EBGK) collision. 
Solution of (|2.6p must be found at every time step and lattice site. Entropic equilibria (also derived from 
the f/-theorem) are always used for ELBGK. 



In the continuous case the Boltzmann H the Maxwellian distribution maximizes entropy and therefore also 
has zero entropy production. In the context of lattice Boltzmann methods a discrete form of the 77-theorem 
has been suggested as a way to introduce thermodynamic control to the system [5] . 

From this perspective the goal is to find an equilibrium state equivalent to the Maxwellian in the contin- 
uum which will similarly maximize entropy. Before the equilibrium can be found an appropriate H function 
must be known for a given lattice. These functions have been constructed in a lattice dependent fashion in 
[8], and H = —S with S from (|2.5|) is an example of a H function constructed in this way. 

Using equilbria derived from a H function with entropy considerations in mind leads to a thermodynam- 
ically correct LBM. This is easy to see in the case of the EBGK collision operator (|2.4p with explicit local 
equilibrium. EBGK collision obvioulsly respect the Second Law (if /3<1), and simple analysis of entropy 
dissipation gives the proper evaluation of viscosity. 

ELBGK finds the value of a that with (3= 1 (inviscid fluid) would give zero entropy production, therefore 
making the position of zero entropy production the limit of any relaxation. For the fixed a used in the 
LBGK method it remains possible, particularly for low viscosity fluids, to relax past this point resulting in 
negative entropy production, violating the Second Law. 

Near to the zero-viscosity limit the LBGK method produces spurious oscillations around Shockwaves. 
Apart from the thermodynamic benefits of using ELBGK it has been claimed [1] that ELBGK's thermo- 
dynamic considerations act as a regularizer. This claim seems to be at odds with other numerical methods 
which respect the same thermodynamic laws as ELBGK. For example the results of Tadmoor and Zhong 
|13) for an entropy correct method display intensive post-shock oscillations. Furthermore it has been demon- 
strated [TOlfTTj that such dispersive oscillations arc artifacts of the lattice rather than thermodynamic issues. 
As ELBGK clearly operates on exactly the same lattice as LBGK and other finite difference schemes it 
warrants a deeper investigation into exactly how it achieves the regularization properties claimed. 

4 Computation of entropic involution 

In order to investigate the stabilization properties of ELBGK it is necessary to craft a numerical method 
capable of finding the non-trivial root in (|2.6j) . In this section we fix the population vectors f and f*, and 




(2.5) 



S(f) = S(f+a(f*-f)). 



(2.6) 



3 The //-theorem for LBMs 



3 



are concerned only with this root finding algorithm. We recast (|2.6[) as a function of a only: 

F(a) = S({ + a{{*-{))-S(t). (4.7) 

In this setting we attempt to find the non-trivial root r of (|4.7[) such that F(r) = 0. It should be noted that 
as we search for r numerically wc should always take care that the approximation we use is less than r itself. 
An upper approximation could result in negative entropy production. 

The following theorem gives cubic convergence order for a simple algorithm for finding the roots of a 
concave function based on local quadratic approximations to the target function. Analogously to the case 
for Newton iteration, the constant in the estimate is the ratio of third and first derivatives in the interval of 
iteration. 

Theorem. For a three times continously differentiable concave entropy function F(a) an iterative root 
finding method based on the zeros of a second order Taylor parabola has cubic convergence sufficiently close 
to the root. 

Proof: Assume that we arc operating in a neighbourhood r G in which F' is negative (as well of course 
F" is negative). At each iteration the new estimate for r is the greater root of the parabola P, the second 
order Taylor polynomial at the current estimate, 

P(a) = F(a n ) + (a- a n )F'(a n ) + (a~ a n f F (4.8) 
The Lagrange remainder form of the error is 

F(a) = F(a n ) + (a-a n )F'(a n ) + {a-a n ) 2 ^^ + (a-a n ) 3 ^p^ 

2 u 

6 

where 7„ lies between a„ and a. Evaluating this at r we see that 

|PW|<i|«„-r| 3 sup|F"'(a)|. 

Now, using the mean value theorem, for some value b n € [on+ijr], 

\P(r)\ = \P(a n+1 ) + (r-a n+1 )P\b n )\>\(r-a n+1 )\M \(F'(b))\. 

Combining the last two equations we see that 

|(r-a n+ i)| < C\a n -rf, where C =\ sup \F"'(a)\ I inf \F'(b)\ . ■ 

OaeJV / b£N 

We use a Newton step to estimate the accuracy of the method at each iteration: 

\a n -r\™\F(a n )/F'(a n )\. (4.9) 

In fact we use a convergence criteria based not solely on a but on a||f* — f||, this has the intuitive appeal 
that in the case where the populations are close to the local equilibrium AS* = S(f *) — S({) will be small and 
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a very precise estimate of a is unnecessary. We have some freedom in the choice of the norm used and we 
select between the standard L\ norm and the entropic norm. The entropic norm is defined as 

\\r-i\\ f , = -(({* -i),D 2 s\ f Af*-{)), 

where D 2 S\^ is the second differential of entropy at point f*, and (x,y) is the standard scalar product. 

The final root finding algorithm then is beginning with the LBGK estimate xq = 2 to iterate using the 
roots of successive parabolas. If this first initiation step produces non-positive population, then the positivity 
rule H] could be used (instead of the mirror image we choose the closest value of a which gives non-negative 
value of populations) . The same rcgularization rule might be suggested if there exists no root we are looking 
for. In the tests described below, this situation never arose. 

We stop the method at the point, 

|a„-r|-||f*-f||<e. (4.10) 

To ensure that we use an estimate that is less than the root, at the point where the method has converged 
we check the sign of F(a n ). If F(a n ) >0 then we have achieved a lower estimate, if F(a n ) <0 we correct 
the estimate to the other side of the root with a double length Newton step, 

an = an ~ 2 F^y 

At each time step before we begin root finding we eliminate all sites with AS < 10 -15 . For these sites 
we make a simple LBGK step. At such sites we find that round off error in the calculation of F by solution 
of equation can result in the root of the parabola becoming imaginary. We note that in such cases a 
mirror image given by LBGK is effectively indistinct from the exact ELBGK collision. 

We now experimentally study the convergence of the method. The convergence of the bisection method 
is presented for control. For the bisection method we calculate an initial estimate using the root of the 
parabola (|4.8|) with — Whichever side of the root this estimate is on, an estimate for the opposite 
side can be found using a double length Newton step. We then have both an upper and lower estimate for 
the root as required for the beginning of the bisection method. For this test e is set to 10~ 75 . This is the 
maximum accuracy following the bound on AS* of 10~ 15 due to the quadratic nature of F. 

For shock tube test using 800 lattice sites at the 400th iteration step (sec detailed description in Scction[5|), 
Fig. [U shows that the parabola based method required two iterations, but not more than two, at some points 
in a vicinity of shock. In other areas one iteration is sufficient for the desired accuracy. Across the whole 
lattice the entropic norm stipulates a slightly greater number of iterations in both methods. 



5 Shock tube tests 

A standard experiment for the testing of LBMs is the one-dimensional shock tube problem. The lattice 
velocities used arc v= (—1,0,1), so that space shifts of the velocities give lattice sites separated by the unit 
distance. 800 lattice sites are used and are initialized with the density distribution 



p(x) 



1, l<x<400, 
0.5, 401 < x < 800. 



Initially all velocities are set to zero. We compare the ELBGK equipped with the parabola based root finding 
algorithm using the entropic norm with the standard LBGK method using both standard polynomial and 
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Figure 1: Iterations required for convergence to e = 10~ 75 under (|4.10[) using (a) Parabola method with L\ 
norm; (b) Parabola method with entropic norm; (c) Bisection method with L\ norm; (d) Bisection method 
with entropic norm. 
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Figure 2: Density profile of the simulation of the shock tube problem following 400 time steps using (a) 
LBGK with polynomial equilibria [v= (1/3) • 10 -1 ]; (b) LBGK with entropic equilibria \v = (1/3) • 10 -1 ]; 
(c) ELBGK [u={l/3)-lQ- 1 ]; (d) LBGK with polynomial equilibria [^ = 10~ 9 ]; (e) LBGK with entropic 
equilibria [v=lQ- g ]; (f) ELBGK [^ = 10~ 9 ]. 



entropic equilibria. The polynomial equilibria arc given in [T^] : 

/* = \ (1 -Zu + Zv?) , fS = f(l-^),/H I 1 + 3 « + 3u2 ) ■ 

The entropic equilibria also used by the ELBGK are available explicitly as the maximum of the entropy 
function ((23)) . 

/*=^(-3u-l + 2Vl + 3w 2 ), / * = ^(2-v / T+3^), /; = ^(37.-1 + 2^1 + 3^). 
Now following (|2.3|) the governing equations for the simulation are 

f-{x,t+l) = f-(x + l,t) + a0(fl(x+l,t)-f-(x + l,t)), 
f (x,t+l) = f (x,t) + aP(fZ(x,t)- f (x,t)), 
U(x,t+l) = f + (x-l,t) + ap(fl(x-l,t)-f + (x-l,t)). 

From this experiment we observe no benefit in terms of regularization in using the ELBGK rather than the 
standard LBGK method (Fig. [2]). In both the medium and low viscosity regimes ELBGK fails to supress 
the spurious oscillations found in the standard LBGK method. 
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Figure 3: Density profile of the simulation of the shock tube problem following 400 time steps using (a) 
LBGK with entropic equilibria and one point median filtering [v= (1/3) • 10 -1 ]; (b) LBGK with entropic 
equilibria and one point median filtering [^ = 10~ 9 ]. 



To explain previous results showing regularization by the ELBGK we note that in the collision integral 
(|2.4p that a and (3 are composite. In this sense entropy production controlled by a and viscosity controlled 
by (3 are the same thing. A weak lower approximation to a would lead effectively to addition of entropy 
at the mostly far from equilibrium sites and therefore would locally increase viscosity. This numerical 
viscosity could, probably, explain the regularization and smoothing of the shock profile seen in some ELBGK 
simulations. 



6 One-Point Median Filtering 

Finally we consider regularizing the LBGK method using median filtering at a single point. We follow the 
prescription detailed in [5]. First, at each time step, we locate the single lattice site x with the maximum 
value of AS(x), and call this value AS X - Secondly, we find the median value of A5* in the three nearest 
neighbours of x including itself, calling this value AS me d- Now instead of being updated using the standard 
BGK over-relaxation this single site is updated as follows: 



f_( x ,t+i)=fi(x+i,t)+J-^{f-(x+i,t)-r_(x+i,t)), 



Mx,t + l) = fS(x,t) + yj-^-{f {x,t)-fS(x,t)), 

f + (x,t+i)=r + (x-i,t)+^^^(u(x-i > t)-fi(x-i > t)). 

We observe that filtering a single point at each time step still results in a significant amount of regularization 

(Fig. ED. 

We also examine in each case the lattice site where the filtering is applied. The zero position is defined 
as the rightmost lattice site with AS > at each time step and the position of the filtering is measured 
relative to this site. The occurrences at each relative position are then summed over the experiment. Wc 
can see (Fig. [4} that the majority of filtering takes place on the shock. However, in the low viscosity case, 
we observe that at a small number of time steps the filtered site moves significantly 'behind' the shockwavc. 
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Figure 4: Distribution of median filtering sites relative to the position of the shock following 400 time steps 
using (a) LBGK with entropic equilibria and one point median filtering ]y= (1/3) • 10 -1 ]; (b) LBGK with 
entropic equilibria and one point median filtering \v= 10~ 9 ]. 

7 Conclusion 

We present three main conclusions from this study. 

1. We do not find any evidence that maintaining proper balance of entropy regularize spurious oscillations 
the Lattice Boltzmann method. For ELBGK we confirm the conclusions of Lax [TU] and Levermore 
with Liu that dispersive oscillations are unavoidable in numerical simulation of shocks. 

2. In order to clean up the parasite dispersive oscillations in the Lattice Boltzmann method it is nec- 
essary to filter the entropy in some way, so as to reduce the extremely-localised incidents of high 
non-equilibrium entropy; see [5]- Previously reported smoothing of shocks must have been via the 
inadvertent introduction of numerical dissipation. (Perhaps, this conclusion could be extended to all 
known regulariscrs of LBM, including those proposed by ourselves in [3].) 

3. For the ID shock tube, one only needs to filter the entropy at one point per time step (usually very local 
to the shock), even at very low viscosity, in order to effectively eliminate the post-shock oscillation. 
We can expect that in 2D and 3D shocks it will be also necessary to filter nonequilibrium entropy in 
some local maxima points near the shock front only. The entropy filtering for non-entropic equilibria 
is possible [5] with use of the Kullback-Leibler distance from current distribution to equilibrium (the 
relative entropy). 

The Matlab code used to produce these results is provided in the appendix. 
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A Matlab Code 



%LBM Lattice Boltzmann Method 

% LBM (MC, V, NX, TM, MOV, ESPIL, LIMIT, NORM) solves the shock tube problem 

% with number of lattice 

% sites NX with separation of one and viscosity V upto TM time steps of 

% length one using method 

% MC : 

% — LBGK polynomial equilibria, 

% 1 — LBGK entropic equilibria, 

% 2 — ELBM Newton iterations, 

% 3 — ELBM Bisection method 

% and entropic limiter 

% LIMIT: 

% - No limiter 

% 1 — Median filtering 

% Output of a movie of the simulation can be controlled with 

% MOV: 

% — No movie, 

% 1 — With movie 

% The accuracy of the root finding for ELBM can be controlled with 

% EPSIL in all cases the root found will result in entropy production. 

% Convergence is partly based on | | feq — f||, the choice of the norm is 

% controlled with 

% NORM: 

% — LI norm 

% 1 — Entropic Norm; 

function [rho, mov, convergence, relativeFiltering, alpha] = ... 

LBM (methodChoice, viscosity, nx, timeMax, MOV, ELBMEpsilon, Limiter, Norm) 

% Over— relaxation parameter 
beta = 1/ (2*viscosity+l) ; 
time = 1; %Start Time 

rho = LBMInitialise (nx) ; %Densities at each lattice site are intialised as 
%the standard shock tube problem 

u = 0; %Velocities at each lattice site are initialised as zero 

populations = LBMQuasiEquilibria (rho, u, methodChoice, nx) ; %Left moving, 
%central and right moving populations at each lattice site respectively are 
%initialised as the quasiequilibria 

relativeFiltering = zeros ( 1 , 25 1 ) ; 

if (MOV == 1) 

mov = moviein (timeMax) ; 
x = 1 : nx; 

else 

mov = ; 

end 

while time < timeMax 

populations = LBMPropagate (populations, nx) ; %Propagate the populations 
%keyboard 



n 



[rho, u] = LBMLatticeParameters (populations) ; %Density and velocity are 
% calculated using the populations 

popequilibriums = LBMQuasiEquilibria (rho, u, methodChoice, nx) ; %New 
% quasiequilibria are found using new density and velocity values 

[alpha , convergence] = LBMEntropicParameter (populations, ... 
popequilibriums, methodChoice, nx, ELBMEpsilon, Norm) ; 

% Finding the non trivial root for constant entropy in the ELBM, for 
% normal LBGK alpha = 

[ limiterSites, shockRelative] = LBMLimiterSites (populations, . . . 
popequilibriums, nx, Limiter) ; 

relat iveFiltering (250 — shockRelative) = ... 

relativeFiltering ( 250 — shockRelative) + 1; 

populations = LBMCollide (populations, popequilibriums, beta, alpha, nx, . . . 
Limiter, limiterSites) ; 

% Populations are collided with the quasiequilibria 

if (-"isreal (populations) ) 
keyboard 
return 

end 

if (MOV == 1) % If movie parameter is enabled record a frame 
unlimitedSites = setdif f ( 1 : nx, limiterSites ) ; 

plot (unlimitedSites, rho (unlimitedSites) , ' . 1 , limiterSites, . . . 

rho (limiterSites) , ' r* ' ) ; 
axis ( [0 nx 0.3 1.3]) 
mov ( : , time) = getframe; 

end 

time = time + 1; ^Increment time 

end 



%A function to intialise the lattice densities for the shock tube problem 

function rho = LBMInitialise (nx) 

rho = zeros ( 1 , nx) ; 

half = floor (nx/2) ; 

rho(l:half) = 1 ; 

rho (half+1 :nx) = 0.5; 



%A function to compute given a current density vector either the polynomial 
%or entropic quasiequilbria . 

function equilibria = LBMQuasiEquilibria (rho, u, methodChoice, nx) 
equilibria = zeros (3, nx); 

if ( methodChoice == 0) %Polynomial Equilibria 

equilibria ( 1 , :) = rho ./6.*(1 — 3*u + 3*u . " 2 ) ; 
equilibria(2,:) = 2. *rho . /3 . * (1 — 3*u."2./2); 
equilibria(3,:) = rho . / 6 . * ( 1 + 3*u + 3*u . " 2 ) ; 
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else %Entropic Equilibria 

equilibria (1, :) = rho . / 6 . * { — 3 . *u — 1 + 2 *sqrt (1 + 3 . *u . " 2 ) ) ; 

equilibria (2,:) = 2. *rho . /3 . * (2 — sqrt (1 + 3 . *u . ~ 2 ) ) ; 

equilibria (3, : ) = rho./6.*( 3 . *u — 1 + 2*sqrt(l + 3.*u."2)); 

end 



%A function to propagate the populations 

function populations = LBMPropagate {old/populations, nx) 

populations = oldpopulations ; 

left = populations (1,1); 

right = populations ( 3 , nx) ; 

populations (1,1: nx— 1) = populations (1,2: nx) ; 
populations (3,2: nx) = populations (3,1: nx— 1 ) ; 
populations (3, 1) = left; 
populations ( 1 , nx) = right; 



%A function to find density and velocity at each lattice site from the 
%populat ions 

function [rho, u] = LBMLatticeParameters (populations) 

rho = populations (1, : ) + populations (2, : ) + populations (3, : ) ; 

u = (populations (3, : ) — populations (1, : ) ) . /rho; 



% Function to find the constant entropy parameter for ELBM, in the case of 
% LBKG this is simply 2 

function [alpha, convergence] = LBMEntropicParameter (populations, ... 
popequilibriums, methodChoice, nx, ELBMEpsilon, Norm) 

%Choice of method: 
%0 — LBGK polynomial equilibria, 
%1 — LBGK entropic equilibria, 
%2 — ELBM Parabola iterations, 
%3 — ELBM Bisection method 

if (methodChoice == ) 

alpha = 2*ones (1, nx) ; 

convergence = zeros (l,nx); 
elseif (methodChoice == 1) 

alpha = 2*ones (1, nx) ; 

convergence = zeros (l,nx); 
elseif (methodChoice == 2) 

convergence = zeros (1, nx) ; 

alpha = 2*ones (1, nx) ; 

Sf = LBMEvaluateS (populations, popequilibriums, zeros ( 1 , nx) ) ; 

SEquil = LBMEvaluateS (popequilibriums, popequilibriums, zeros ( 1 , nx) ) ; 

remaining = find (abs ( SEquil — Sf) > 10"— 15); 

popnorm (remaining) = LBMNorms (populations ( : , remaining) , ... 
popequilibriums ( : , remaining) , length (remaining) , Norm) ; 

iteration = 1; 

while (isempty (remaining) == 0) 
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alpha (remaining) = LBMInteriorParabola (populations ( : , remaining) , 
popequilibriums ( : , remaining) , Sf (remaining) , alpha (remaining) ) 

Spop = LBMEvaluateS (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , alpha (remaining) ) ; 
Sdif f = LBMEvaluateDif f S (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , alpha (remaining) ) ; 

AAlpha ( remaining) = (Spop — Sf (remaining) ) ./Sdiff; 

nowdone = find( popnorm (remaining) . *abs (AAlpha ( remaining) ) ... 

< ELBMEpsilon) ; 
nowremaining = f ind ( popnorm ( remaining) .* ... 

abs (AAlpha (remaining) ) > ELBMEpsilon) ; 

done = remaining (nowdone) ; 
remaining = setdif f (remaining, done) ; 

SDone = LBMEvaluateS (populations ( : , done) , ... 

popequilibriums ( : , done) , alpha (done) ) — Sf (done) ; 
negFind = findfSDone < 0); 
negativeEntropy = done (negFind) ; 

alpha (negativeEntropy) = alpha (negativeEntropy) . . . 
— 2*AAlpha (negativeEntropy) ; 

SDone = LBMEvaluateS (populations ( : , done) , ... 

popequilibriums ( : , done) , alpha (done) ) — Sf (done) ; 
if (SDone < 0) 

keyboard 

end 

convergence (done) = iteration; 

iteration = iteration + 1; 
if (-dsreal (alpha) ) 
keyboard 

end 

if (iteration > 200) 
keyboard 

end 

end 

elseif (methodChoice == 3 ) 

convergence = zeros (l,nx) ; 
alpha = 2*ones (1, nx) ; 

Sf = LBMEvaluateS (populations, popequilibriums, zeros ( 1 , nx) ) ; 

SEquil = LBMEvaluateS (popequilibriums, popequilibriums, zeros ( 1 , nx) ) ; 

remaining = find (abs ( SEquil - Sf) > 10"— 14); 

popnorm (remaining) = LBMNorms (populations ( : , remaining) , ... 
popequilibriums ( : , remaining) , length (remaining) , Norm) ; 

alpha (remaining) = LBMInteriorParabola (populations ( : , remaining) , . . . 

popequilibriums ( : , remaining) , Sf (remaining) , alpha (remaining) ) ; 
Spop = LBMEvaluateS (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , alpha (remaining) ) — Sf (remaining) ; 
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leftside = find(Spop > 0); 
rightSide = find(Spop < 0); 

left (remaining (leftside) ) = alpha (remaining (leftside) ) ; 
right (remaining (rightSide) ) = alpha (remaining (rightSide) ) ; 

Spop = LBMEvaluateS (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , alpha ( remaining) ) ; 
Sdif f = LBMEvaluateDif f S (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , alpha (remaining) ) ; 

AAlpha (remaining) = {Spop — Sf (remaining) ). /Sdif f; 

right (remaining (leftside) ) = alpha (remaining (leftside) ) . . . 

— 2*AAlpha (remaining (leftside) ) ; 
left (remaining (rightSide) ) = alpha (remaining (rightSide) ) . . 

+ 2*AAlpha (remaining (rightSide) ) ; 

if (right < left) 
keyboard 

end 

iteration = 1 ; 

mid = left + (r ight-lef t ) 12 ; 



while (isempty (remaining) == 0) 

mid ( remaining) = left (remaining) . . . 

+ (right (remaining) —left (remaining) ) /2; 



Sleft = S (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , left (remaining) ) ; 
Sright = LBMEvaluateS (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , right (remaining) ) ; 
Smid = LBMEvaluateS (populations ( : , remaining) , ... 

popequilibriums ( : , remaining) , mid (remaining) ) ; 

nowdone = find( popnorm (remaining) . * (right (remaining) ... 

— left (remaining) ) < ELBMEpsilon) ; 

nowremaining = f ind ( popnorm ( remaining) . * (right (remainin 1 

— left (remaining) ) > ELBMEpsilon); 
done = remaining (nowdone) ; 
remaining = setdif f (remaining, done) ; 
alpha (done) = left (done) ; 
convergence (done) = iteration; 

for itr = 1 : length ( remaining) 

if (Smid (nowremaining (itr) ) > Sf (remaining (itr) ) ) 
left (remaining (itr) ) = mid (remaining (itr) ) ; 

else 

right (remaining (itr) ) = mid (remaining (itr) ) ; 

end 

end 

iteration = iteration + 1; 
if (-"isreal (alpha) ) 
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keyboard 

end 

if (iteration > 100) 

keyboard 

end 

end 

end 



% A function to evaluate the entropy. 

function S = LBMEvaluateS (populations, popequilibriums, alpha) 
alphapop = populations + (ones ( 3 , 1 ) *alpha) . * (popequilibr iums ... 

— populations) ; 

S = —alphapop ( 1 , : ) . *log (alphapop ( 1 , : ) ) ... 

— alphapop ( 2 , : ) . *log (alphapop (2,:)./4) ... 

— alphapop (3, : ) . *log (alphapop (3, : ) ) ; 



% A function to implement the norms necessary to measure convergence of the 
% root finding 

function norms = LBMNorms (populations, popequilibriums, nx, nChoice) ; 
norms = zeros (l,nx); 
if (nChoice == 0) 
for j = 1 : nx 

norms (j) = norm (populations (:, j ) — popequilibriums ( : , j ) , 1 ) ; 

end 

elseif (nChoice == 1) 
for j = 1 : nx 

poptemp = (populations ( 1 , j ) ... 

— popequilibriums ( 1 , j ) ) . " 2 . /popequilibriums ( 1 , j ) ; 
poptemp = poptemp + (populations ( 2 , j ) ... 

— popequilibriums ( 2 , j ) ) . " 2 . /popequilibriums ( 2 , j ) ; 
poptemp = poptemp + (populations ( 3 , j ) ... 

— popequilibriums (3, j) ) . "2 . /popequilibriums (3, j) ; 
norms (j) = sqrt (poptemp) ; 

end 

end 



% A function to find an interior approximation to the root of the entropy 
% parabola 

function intPab = LBMInteriorParabola (populations, popequilibriums, ... 
STarget, alpha) 

SAlphaZero = LBMEvaluateS (populations, popequilibriums, alpha) ; 
SDashAlphaZero = LBMEvaluateDif f S (populations, popequilibriums, alpha) ; 
SDash2AlphaZero = LBMEvaluateDif f2S (populations, popequilibriums, alpha) ; 
intPab = [ ] ; 
for j = 1 : length (alpha) 

intPab(j) = max (roots ( [ . 5 . *SDash2AlphaZero ( j ) SDashAlphaZero ( j ) ... 
SAlphaZero(j) — STarget ( j )]) ) + alpha(j); 

end 



% A function to find the first derivative of the entropy 
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function Dif f S = LBMEvaluateDif f S (populations, pop equilibriums, alpha) 
alphapop = populations + (ones (3,1) *alpha) . * (popequilibriums . . . 

— populations) ; 
partDiff = popequilibriums — populations; 
Dif f S = -partDiff (1, : ) . * (log (alphapop (1, : ) ) +1) ... 

-partDiff (2, :).* (log (alphapop (2, :). /4) + 1) ... 

-partDiff (3, :).* (log (alphapop (3, :) ) + 1 ) ; 



% A function to find the second derivative of the entropy. 

function Dif f 2S = LBMEvaluateDif f 2S (populations, popequilibriums, alpha) 

alphapop = populations + (ones (3, 1) *alpha) . * (popequilibriums — populations); 

partDiff = popequilibriums — populations; 

Dif f 2S - - partDiff (1, : ) . "2 . / (alphapop (1, : ) ) ... 

- partDiff (2, : ) . "2 . / (alphapop (2, : ) ) ... 

- partDiff (3, : ) . "2 . / (alphapop (3, : ) ) ; 



% A function to detect an appropriate lattice site for limiting and to 
% measure it's position relative to the leading edge of the shock, 
function [Sites, shockRelative] = LBMLimiterSites (populations, . . . 

popequilibriums, nx, Limiter) 
shockRelative = 0; 
if (Limiter == 0) 

Sites = [ ] ; 
elseif (Limiter == 1} 

Sf = LBMEvaluateS (populations, popequilibriums, zeros ( 1 , nx) ) ; 

Sfeq = LBMEvaluateS (popequilibriums, popequilibriums, zeros (1, nx) ) ; 

Deltas = Sfeq - Sf; 

shockPos = max(find(DeltaS > 1CT-15) ) ; 
[M, I] = max (DeltaS) ; 
Sites = I; 

shockRelative = shockPos — Sites; 

end 



%Function for simple BGK collision of populations with quasiequilibria with 
% a limiter applied if necessary 

function populations = LBMCollide (oldpopulations, popequilibriums, beta, . . . 

alpha, nx, Limiter, limiter Sites) 
unlimitedSites = setdif f ( 1 : nx, limiter Sit es ) ; 

populations ( : , unlimitedSites) = oldpopulations ( : , unlimitedSites) . . . 
+ beta . * (ones (3,1) * alpha (unlimitedSites) ) ... 
. * ( popequilibriums ( : , unlimitedSites) . . . 
— oldpopulations ( : , unlimitedSites) ) ; 
if (Limiter == 1) 

Sf = LBMEvaluateS (oldpopulations (:, (limiterSites — 1): ... 

(limiter Sites + 1 ) ) , popequilibriums ( : , (limiterSites — 1 ) : ... 

(limiterSites + 1) ) , zeros(l, 3) ) ; 
Sfeq = LBMEvaluateS (popequilibriums ( : , (limiterSites — 1) : ... 

(limiterSites + 1 ) ) , popequilibriums ( : , (limiterSites — 1 ) : ... 

(limiterSites + 1) ) , zeros (1, 3) ) ; 
AS = Sfeq - Sf; 
Smed = median (aS); 
coeff = sqrt (Smed/AS (2) ) ; 
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populations ( : , limiterSites) = popequilibriums ( : , limiterSites ) 
+ coef f . * (oldpopulations ( : , limiterSites) . . . 
— popequilibriums ( : , limiterSites) ) ; 

end 
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